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We present a general numerical approach to construct local Kohn-Sham potentials from 
orbital-dependent functionals within the all-electron full-potential linearized augmented-plane-wave 
(FLAPW) method, in which core and valence electrons are treated on an equal footing. As a prac- 
tical example, we present a treatment of the orbital-dependent exact-exchange (EXX) energy and 
potential. A formulation in terms of a mixed product basis, which is constructed from products of 
LAPW basis functions, enables a solution of the optimized-effective-potential (OEP) equation with 
standard numerical algebraic tools and without shape approximations for the resulting potential. 
We find that the mixed product and LAPW basis sets must be properly balanced to obtain smooth 
and converged EXX potentials without spurious oscillations. The construction and convergence of 
the exchange potential is analyzed in detail for diamond. Our all-electron results for C, Si, SiC, 
Ge, GaAs semiconductors as well as Ne and Ar noble-gas solids are in very favorable agreement 
with plane-wave pseudopotential calculations. This confirms the adequacy of the pseudopotential 
approximation in the context of the EXX-OEP formalism and clarifies a previous contradiction 
between FLAPW and pseudopotential results. 
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I. INTRODUCTION 

Its wide applicability, accuracy, and computational ef- 
ficiency have made density-functional theory (DFT)^'^ 
the standard method for describing the ground state of 
many-electron systems. The vast majority of practical 
calculations employ the Kohn-Sham (KS) formalism;^ 
where the interacting many-electron system is mapped 
onto an auxiliary noninteracting system. In this KS sys- 
tem the noninteracting electrons move in a local effective 
potential that is defined in such a way that the electron 
densities of the real and auxiliary systems coincide. This 
potential is the sum of the external, Hartree, and the 
exchange-correlation (xc) potentials. The latter two con- 
tributions take into account the electron-electron interac- 
tions including, in an indirect way, all many-body effects. 
The form of the density-functional for the xc potential, 
which can further be divided into an exchange and a cor- 
relation term, is unknown and must be approximated in 
practice. 

Fortunately, already simple approximations, like the 
local-density approximation (LDA)^iS, or generalized gra- 
dient approximation (GGA),^'^ give reliable results for a 
wide range of materials and properties. Nevertheless, the 
LDA and GGA suffer from several shortcomings. First, 
the electrostatic interaction of the electron with the total 
electron charge, described by the Hartree potential, con- 
tains an unphysical interaction of the electron with itself, 
commonly referred to as Coulomb self-interaction. This 
extra term should be compensated exactly by an iden- 
tical term with opposite sign in the exchange potential, 
in the same manner as in Hartree-Fock theory. However, 



as the LDA and GGA exchange potentials are only ap- 
proximate, this cancellation is incomplete and part of the 
self-interaction remains. This error leads, in particular, 
to an improper description of localized states, which ap- 
pear too high in energy and tend to delocalize. Second, 
the LDA and GGA functionals do not give rise to a dis- 
continuity of the xc potential with respect to changes in 
the particle number. This discontinuity should in general 
be finite (and positive), as it corresponds to the differ- 
ence of the real and the KS band gap.^'^ The latter is 
well known to underestimate experimental gaps by typi- 
cally 50% or even more. This is often called the band-gap 
problem of LDA and GGA. The significance of the dis- 
continuity for a meaningful prediction of the fundamental 
band gap is discussed in Refs. 1(1 and [nl. 

Functionals that depend explicitly on the electron or- 
bitals and thus only implicitly on the electron density 
form a new generation of xc functionals ji^i^ Already the 
simplest variant, the EXX functional)^"— which treats 
electron exchange exactly but neglects correlation alto- 
gether, remedies the aforementioned shortcomings of the 
more conventional local and semilocal functionals: the 
Coulomb self-interaction term is exactly canceled, and 
the local EXX potential exhibits a nonzero discontinu- 
ity at integral particle numbers because of the orbital 
dependence. 

After applications to atomsii"— the first implemen- 
tation for periodic systems was published in 1994 by 
Kotani^*^ who employed the atomic sphere approxima- 
tion within the linearized muffin-tin orbital method. In 
this approximation only the spherical part of the poten- 
tial around each atom is taken into account. The KS 
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band gap turned out to be closer to experiment than in 
LDA or GGA.2°'2^ This indicates that the xc disconti- 
nuity and the effect of neglecting electron correlation for 
these systems are roughly of the same magnitude but of 
different sign and thus tend to cancel each other. Later, 
results even closer to experimental values were obtained 
from plane-wave calculations—^— employing the pseu- 
dopotential (PP) approximation, which allows accurate 
treatment of the warped shape of the EXX potential ex- 
cept for the regions close to the atomic nuclei where the 
potential is smoothed. 

However, when the first all-electron (AE), full- 
potential results were reported, they deviated substan- 
tially from the PP values and were also in considerably 
worse agreement with experiment. In Ref. [25|, Sharma 
et al, who implemented the EXX functional within the 
FLAPW method, argue that the success of EXX in the 
earlier PP calculations is only an artifact of the neglect 
of the core-valence exchange interaction. They conclude 
that treatment of core and valence electrons on the same 
footing is imperative for a proper EXX calculation. This 
work started a controversy about the adequacy of the 
PP approximation in the EXX formalism. Recently, 
Engel^^ showed that, on the contrary, AE and PP re- 
sults for lithium and diamond differ only marginally. The 
AE calculations were performed with a plane-wave basis 
by pushing the PP and plane-wave cutoffs to the AE 
limit. Also recently, Makmal et al. reported a similarly 
good agreement between AE and PP calculations for the 
diatomic molecules BeO and CO using real-space grid 
approaches i21 

A final comparison between PP results and results ob- 
tained from a genuine AE approach for periodic systems, 
such as the FLAPW method, is still missing, though. 
With this work we want to fill this gap. We present an 
alternative implementation of the EXX approach within 
the FLAPW method, which uses a numerical approach 
different from the one reported in Ref. [2^ It employs a 
specifically designed basis, the mixed product basis, in 
which the optimized-efTective-potential (OEP) equation 
for the EXX potential is solved. The mixed product func- 
tions form an all-electron basis for the products of single- 
particle wave functions occurring in the OEP equation. 
Within this approach, both spherical and nonspherical 
as well as warped interstitial contributions to the EXX 
potential are fully taken into account. For the exam- 
ple of diamond, we discuss in detail the convergence of 
the local EXX potential and the resulting KS band gaps 
with respect to the quality of the mixed product and 
LAPW basis sets. We demonstrate that a smooth po- 
tential and a direct KS band gap very close to the result 
by Engel^*' are obtained if the two basis sets are prop- 
erly balanced. A similar behavior has been reported for 
Gaussian and plane-wave basis setsi^"— To validate our 
findings, we also report results for Si, SiC, Ge, GaAs, 
and crystalline Ne and Ar. For all materials, we find 
a very good agreement between our AE and previously 
published plane- wave PP values. We conclude that the 



large discrepancies found in Ref. |25| cannot be attributed 
to the core-valence interaction. 

The paper is organized as follows. Sections [IT] and IIIII 
give brief introductions into the theory and the FLAPW 
method. Our implementation of the EXX functional 
within the FLAPW program Fleur'^^ is described in 
Section IIVI Section |V] discusses the convergence of the 
effective potential for the example of diamond and com- 
pares AE EXX KS eigenvalue differences commonly in- 
terpreted as transition energies for C, Si, SiC, Ge, GaAs, 
crystalline Ne, and Ar with theoretical plane-wave and 
experimental values from the literature. Finally, we draw 
our conclusions in Section [VIl 



II. THEORY 



The KS formalism"^ of DFT^ relies on an auxiliary sys- 
tem of noninteracting electrons, which move in the spin- 
dependent local effective potential 



V,Ur) = VMv) + Vnir) + VZir) 



(1) 



with the external, Hartree, and xc potential, respectively. 
The latter is defined in such a way that the electron spin 
densities coincide with those of the real interacting sys- 
tem. It is given by the functional derivative of the xc 



energy functional Ey^ 



i^] with respect to the electron 



spin density n'^(r) {a =tii): 
K'^c(r) 



(2) 



The orbitals describing the electrons in the auxiliary sys- 
tem obey the KS equations 
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(3) 



where denotes the KS orbital of spin cr, band index n 
and Bloch vector k. Hartree atomic units are used except 
where explicitly noted. The electron spin density is given 
by a sum over the occupied states 



occ. 



(4) 



By a summation over Bloch vectors we mean an integra- 
tion over the Brillouin zone, which is sampled by a finite 
set of mesh points. 

For the conventional LDA and GGA functionals, the 
xc energy functional depends locally on the spin densities 
and, in the case of GGA, on the their gradients, and the 
functional derivative in Eq. Q translates to a derivative 
of a function and is evaluated in a straightforward way. 
However, for orbital-dependent functionals, E^c depends 
only indirectly on the electron spin densities: the KS 
orbitals, which define i?xc [v'^i 'P^] j ^re functionals of the 
effective potential Vj,g(r) through Eq. ([3]), and V^ff{r) is 
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a functional of n" . Therefore, one must apply the chain 
rule to calculate the functional derivative in Eq. ^ 



K'^c(r) 



(5) 



where the sum runs over all KS states present in E^c- 
Then, multiplication with the single-particle spin-density 
response function 



xr(r,r') = 



(5n'"(r) 



(6) 



integration, and use of x^(r, r') — x^(r',r) yields an in- 
tegral equation for the xc potential 



x:{ry)v:Ar')dV 
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(7) 



In this work we employ, as a practical example, the 
orbital-dependent EXX functional 



(7 k,q n.n' 

<k(r)<'q(r)<rq(r')¥>;^k(r') 



(8) 



whose functional derivative with respect to the KS wave 
functions is given by the well-known Hartree-Fock ex- 
pression 



5E^ 



^^k(r")KV(r",r')d-^r" (9) 



with 



q n' ' 

First-order perturbation theory yields the wave-function 
response 



'^<k(i-) 



E 



<'k(r'Xk(r') 



^^'k(r) (11) 



and together with Eq. ^ the spin-density response func- 
tion 



occ. unocc. 



<k(r)¥':^'k(rX'*k(r')<k(r') 
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where time-reversal symmetry has been used. Using 
Eqs. dS]), m]), and the integral equation [Eq. ^\ 
turns into 



with 



r(r) 



X:(r,r')Kr(r')d^r' = ^(r) 



(13) 
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k n 



\^(i-)<k(r) 



^Tik ^n'k 



(14) 



and 



KklKVK'i 



(15) 



/ \ <k(r)KV('^''^')<'k(r')rf''^d3r'. 



In this form the integral equation is called OEP equa- 
tion and goes back to Sharp and Horton,"^^ who derived 
Eq. (I13p as a result of a variational minimization of the 
Hartree-Fock total energy under the additional constraint 
that the orbitals experience a local rather than a nonlo- 
cal potential. Sahni et aL— finally realized that the OEP 
approach is equivalent to the construction of a local EXX 
potential within the KS formalism. 



III. FLAPW METHOD 

The LAPW basis^"— is constructed from piecewise 
defined functions to deal, at the same time, with the 
atomic-like potential close to the nuclei and the smooth 
potential in the region far away from the nuclei. For this 
purpose, space is partitioned into nonoverlapping atom- 
centered muffin-tin (MT) spheres and the remaining in- 
terstitial region (IR), where the smoothness of the po- 
tential allows to employ plane waves as basis functions. 
At the MT sphere boundaries, these plane waves are 
matched in value and first radial derivative to linear com- 
binations of spin-dependent MT solutions u^q {r)Yim{^) 
of the radial scalar-relativistic Dirac equation and their 
energy derivatives uf^ {r)Yim(i) using the spherical aver- 
age of the effective potential and predefined energy pa- 
rameters that lie in the energy range of the occupied 
states. Here, Y/m(f) denote the spherical harmonics, r 
is measured from the MT center of atom a and r = r/r 
is a unit vector. This gives the LAPW basis functions 
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^exp[z(k + G).r] 



ifr e IR 
if r e MT(a) 



(16) 



r 



for the valence electrons with the unit-cell volume 
and reciprocal lattice vectors G. For a practical cal- 
culation cutoff values for the reciprocal lattice vectors 
|k+ G| < Gmax and the angular momentum I < Zmax are 
employed. The core states are obtained by solving the 
fully relativistic Dirac equation with the spherical aver- 
age of the effective potential. 

The basis functions, defined in Eq. (jl6p . can repre- 
sent only those wave functions accurately whose energies 
are sufficiently close to the energy parameters, which are 
usually located in the valence-band region. For a precise 
description of semicore and high-lying unoccupied states 
that are far away from the energy parameters the basis 
must be augmented and local orbitals (lo)^"— are cur- 
rently the best developed technique. Let us assume that 
we want to improve the basis for states with an angular 
momentum I around an energy e'°. Then we construct 
an additional radial function it°'^(r, e'°) from the radial 
scalar-relativistic Dirac equation with the energy pa- 
rameter and form a linear combination uf^ {r)Yimii') 
{p > 2, the index p is a label numbering the basis func- 
tions for a given I, a, and a) from uf'^{r, e'°) and the ra- 
dial functions u'^^ir) and u1'^{r), already defined above, 
such that ufp{r) is normalized and its value and radial 
derivative vanish at the MT boundary. In this way, the 
local orbital u^^ {r)Yim,{r) is completely confined to the 
MT sphere and need not be matched to a plane wave 
outside. For semicore states, which are nearly dispersion- 
less, the energy parameter e'° is fixed at the semicore 
energy level. For the unoccupied states we use energy 
parameters chosen such that the solutions of the radial 
scalar-relativistic Dirac equation fulfill 



(17) 



r=S 



at the MT sphere boundary r = S, following a proce- 
dure proposed in Ref. HO- This condition yields for each 
I quantum number a series of orthogonal solutions of in- 
creasing energies. We use the resulting local orbitals to 
converge the LAPW basis in a systematic way. 



IV. IMPLEMENTATION 

To solve the integral equation [Eq. (fT3|) ]. we introduce a 
basis {M/(r)} that reformulates the equation as a linear- 
algebra problem 



J 



(18) 



which can be solved for the exchange potential by 
matrix inversion of Xs applying standard numerical tech- 
niques. As all quantities appearing in Eq. ()13|) are defined 
in terms of wave-function products, the basis should be 
constructed foremost of products of LAPW basis func- 
tions. In recent publications we have already used such 
a mixed product basis (MPB), which was first proposed 
by Kotani and van Schilfgaarde,**^ to implement hybrid 
functionals^^ and the GW approximation^'^ as well as 
calculate EELS spectra.''^ However, we will introduce 
a slightly modified version for the present purpose: (1) 
since the potential is strictly periodic, the MPB may be 
restricted to k = 0, (2) we add the atomic exact ex- 
change potential as a basis function, and (3) we make 
the functions Mj{r) continuous over the whole space. 

The construction of the MPB and the implementation 
of the spin-density response function Xsij ^'^'^ ^/ 
described in Sec. II V Al and ITVBl respectively. Numerical 
tests of the implementation are shown in Sec. IIV CI 



A. Mixed product bcisis 

The MPB consists of plane waves in the IR and MT 
functions Af£p(r)Y£,M(f ) in the spheres that derive from 
products of the functions uf^ {r)Yim{i)- As in the LAPW 
basis, cutoff values GJnax for the interstitial plane waves 
and Lmax for the angular momentum quantum numbers 
L are employed. For mathematical details of the con- 
struction of the MPB we refer the reader to our previ- 
ous publications Refs. |4^ - |43 . Here, we lay emphasis on 
the modifications for the present implementation of the 
EXX-OEP method. 

From EXX-OEP calculations of atoms it is known that 
the local exchange potential shows pronounced humps 
which reflect the atomic shell structure As the elec- 
tron orbitals contract spatially for atoms with larger 
atomic numbers, these humps move closer and closer to 
the atomic nucleus. Near the nucleus, the exchange po- 
tential of a periodic crystal resembles that of the corre- 
sponding atom, because the long-range exchange inter- 
actions with electrons on neighboring atoms contribute 
only a slowly varying potential there. Therefore, we aug- 
ment the MPB with the spherical exchange potential 
from an atomic EXX-OEP calculation performed with 
the relativistic atomic structure program RELKS4^"— 
It is added to the set of spherical MT functions (L = 0). 
The rest of the basis must then only describe the differ- 
ence between the atomic and the crystal exchange poten- 
tial. In this way all nonlocal exchange contributions are 
fully taken into account. 
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To avoid discontinuities of the resulting potential at 
the MT sphere boundaries, we form linear combinations 
of the MT functions and interstitial plane waves that are 
continuous in value and first derivative there. In anal- 
ogy to the construction of the LAPW basis (s. Sec. IIIip . 
two radial functions per Im channel are used to augment 
the interstitial plane waves in the MT spheres, while the 
remaining functions are combined to form local orbitals. 
We note that there are usually far more than two radial 
functions per Im channel in the MPB. We also note that 
such a construction was not needed in our earlier imple- 
mentations. 



B. Spin-density response function and t/(r) 

In the MPB the spin-density response function in 
Eq. dHD and (r) in Eq. (??) become 



occ. unocc. 



and 



k n 



occ. unocc. 



(19) 



(20) 

Both core and valence states are taken into account in 
the sums over the occupied states in Eqs. (ITO)) and (pOj) . 

As said before, solving Eq. for the exchange po- 
tential involves the matrix inversion of Xs ij- The Ho- 
henberg and Kohn theorem guarantees that the response 
function is invertible except for variations of the poten- 
tial given by an addition of a constant. The latter restric- 
tion gives rise to a constant eigenfunction of XsU with 
eigenvalue 0, which we eliminate from the outset by or- 
thogonalizing all MPB functions to a constant function 
such that variations in the potential by a constant are 
excluded. 

Equation (|20p contains the matrix elements of the 
nonlocal exchange potential [Eq. ([T5|) ] between occupied 
(core and valence) and unoccupied states. In a recent 
publication, we described an efhcient scheme to calcu- 
late the valence-valence and valence-conduction matrix 
elements within the FLAPW method.**^ For the present 
implementation, this scheme has been extended to the 
core-conduction matrix elements. Furthermore, spatial 
and time-reversal symmetries are exploited to restrict the 
k-point sums to the irreducible wedge of the Brillouin 
zone in Eqs. (HH) and (HH])-^'^ 



C. Numerical tests 

In this section, we present numerical tests of the spin- 
density response function, the function i/(r), and the 



Table I: Numerical test for the response function Xs(r,r') — 
Sn{r)/5Vs{r') of diamond. The exact response of the density 
An(r) is compared with its linear approximation An''"(r) = 
J2nTMi{r) with n^r = « E.j XsJJ^pcr,j. The L2 norm 
[/ I An''"(r) — An(r)pd''r]^''^clearly shows a quadratic depen- 
dence on the perturbation strength a. 



0.01 



0.001 



0.0001 



|An""-An| 2.16710 x 10""* 2.17018 x lO"'' 2.17035 x 10^ 



Table II: Same as Table U for t(r) = SE^/SVs{r). The dif- 
ference of the exact response AE^ of the exchange energy 
and its linear approximation A^B^'" — a'^j tiVpcrj depends 



quadratically 


on a. 






Q 


0.01 


0.001 


0.0001 




1.69455 


1.74639 X 10"^ 


1.75171 X 10"^ 


Ai5i-" 


1.75230 


1.75230 X 10"^ 


1.75230 X 10^2 


1A£^-A£,| 


5.77479 X 10-2 


5.90882 X 10"'' 


5.92211 X 10"^ 



resulting exchange potential. According to the Eqs. ([2]), 
([ni, and (??) all three quantities are functional deriva- 
tives of the form j^j^- Thus, they describe the linear 
response of a quantity A with respect to changes of a 
quantity B. 

For the case of diamond, we calculate the changes 
An(r) and A£^x that result from an explicit per- 
turbation Vper(r) = aX]/ Vper,/M/(r), where Vper,/ 
are random numbers, by exact diagonalization of 
the perturbed Hamiltonian and check whether they 
correspond to their linear counterpart An''"(r) = 
fXs{r,r')Vp,,{r')d^r', A^- = f t{r')Vp,,{r')d^r' , and 



A£;^'" = J Vx{r')An{r)d^r up to linear order m 

Tables mini and IIIII show that, indeed, the differences 
I An"" - An I and |A£;^'" - AE'xl depend quadratically on 
the perturbation strength a which confirms the validity 
of our implementation. 



V. RESULTS & DISCUSSION 

In this section, we present results for a variety of semi- 
conductors and insulators obtained with our implemen- 
tation of the EXX-OEP approach within the FLAPW 
method. In particular, we demonstrate for the case of 



Table III: Same as Table HH for 14 (r) 
A£i'" = Ez 14,7 /M;(r)An(r)dV. 



5E^/5n{r) with 



0.01 



0.001 



0.0001 



AE^ 1.69455 1.74639 x 10^^ 1.75171 x 10"^ 

ASi'" 1.67924 1.74479 x 10~^ 1.75155 x 10"^ 

\AEl^-AE^\ 1.53093 x 10"^ 1.59919 x 10"" 1.60625 x 10"® 
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diamond that a smooth and physical EXX potential re- 
quires a balance of the basis sets: the LAPW basis for 
the wave functions must be converged with respect to 
a given MPB until the EXX potential does not change 
anymore. This is somewhat counterintuitive and in con- 
trast to our implementation of the hybrid functionals 
where, conversely, the MPB must be converged for a 
given LAPW basis. A similar behavior has been found in 
implementations employing plane- wave and Gaussian ba- 
sis setsi ^^'^^ We will analyze and explain this point later 
in this section. 

Figure [T] shows the local EXX potential on lines con- 
necting two neighboring carbon atoms along the [111] 
[Fig. [Ha)] and the [100] [Fig.[Ti;b)] directions, see Fig.d 
A 4x4x4 k-point sampling is employed, and the MPB 
parameters are G'^^^ = 3.4 Oq"^ (cq is the Bohr radius) 
and Lmax — 4, giving rise to five s-, four p-, four d-, 
and three /-, and two g-type radial functions per atom. 
These cutoff values are well below those of the LAPW 
basis, ^max = 6 and Gmax = 4.2 Oq^, which reflects the 
relative smoothness of the potential compared with the 
shape of the wave functions. However, if we only use the 
conventional basis of augmented plane waves, defined in 
Eq. (fTB|) . the potential (dashed lines) shows an overpro- 
nounced intershell hump and tends to an unphysical pos- 
itive value close to the atomic nucleus (r ~ 0). This is 
a case where the basis sets are unbalanced. In particu- 
lar, the LAPW basis lacks flexibility in the MT spheres 
as becomes obvious when we add local orbitals, which 
are nonzero only in the MT spheres. We find that six 
local orbitals for each Im channel with I — 0, . . . , 5 and 
|m| < Z, placed at higher energies according to the pre- 
scription described in Sec. IIIIl are needed to converge 
the local EXX potential. This is reasonable since, with 
the cutoff L,„ax = 4, the occupied 2s and 2p states of 
diamond couple maximally to the I = 5 contribution of 
the unoccupied states. With so many local orbitals the 
number of basis functions is increased roughly by a fac- 
tor of five: there are about 100 augmented plane waves 
(the exact number depends on the k point) and 432 addi- 
tional local orbitals. All resulting KS bands, about 530, 
are taken into account in the sums of Eqs. and (fHl) . 
The resulting potential is shown as solid lines in Fig. [1] 
and looks smooth and physical. It is remarkable that, 
even for diamond, it takes so much effort to converge the 
EXX potential since, in conventional LDA or GGA cal- 
culations, diamond is treated readily with a very modest 
LAPW basis without any local orbitals. 

Before analyzing this point in more detail, we want to 
identify the MPB functions that contribute most to the 
MT part of the potential. Not surprisingly, the function 
that corresponds to the atomic EXX potential gets the 
largest weight. In fact, close to the atomic nuclei this 
function (dotted lines in Fig. [T]) and its bulk counterpart 
are indistinguishable. They deviate more and more to- 
wards the MT sphere boundary (5*= 1.42 oq), where the 
atomic EXX potential already enters the typical 1/r be- 
havior, while the crystal EXX potential is periodic. The 




(a) 




(b) 

Figure 1: Local EXX potential of diamond on lines starting 
and ending at atomic nuclei along (a) the [111] and (b) the 
JlOO] directions. Dashed and solid lines correspond to cases 
where the MPB and LAPW basis sets are unbalanced and bal- 
anced, respectively. For comparison, we also show the atomic 
EXX potential of carbon as dotted lines, shifted to align with 
the crystal EXX potential at the atomic nucleus at r = 0. 



second largest contribution comes from the constant MT 
function, which helps to align the MT potential to the 
interstitial one. We note that there is no ambiguity with 
respect to adding a constant to the potential over the 
whole space, since the constant function has been elimi- 
nated explicitly from the MPB (see Sec. lIVEJt giving rise 
to the condition / Vy^{Y)d'^r = 0. 

So far, we have only discussed the MT potential, whose 
proper convergence requires additional local orbitals in 
the spheres. We find an analogous behavior for the inter- 
stitial potential. As is seen in Fig. [31 the cutoff radius of 
the reciprocal lattice vectors included in the LAPW ba- 
sis set must be converged with respect to that of the 
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Figure 2: Cubic unit cell of diamond with the (Oil) plane. 
The lines connecting two carbon atoms along the [111] and 
[100] directions are indicated as solid lines. 
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Figure 3: Convergence of the interstitial EXX potential on a 
line connecting two carbon atoms along the [100] direction for 
four different reciprocal cutoff radii of the LAPW basis. The 
corresponding MPB cutoff is fixed to 5.8 Oq^. 



MPB. To show this effect clearly, the latter was cho- 
sen much larger than necessary, G'^^^ — b.d>aQ^. Sim- 
ilarly to the MT potential, the interstitial potential ex- 
hibits spurious oscillations in the underconverged cases. 
Only if Gmax is large enough, Gmax > 5.0 0(7^ the os- 
cillations are suppressed and a smooth potential is ob- 
tained. Furthermore, Fig. [3] shows that, in the under- 
converged cases, the potential is not continuous at the 
MT sphere boundary because the oscillatory potentials 
possess large- G Fourier coefficients and require spherical 
harmonics beyond L^ax = 4 in the spheres for a proper 
matching. Fortunately, the converged EXX potential is 
a smooth function and already moderate reciprocal cut- 
off radii are sufficient, typically 75% of the usual LAPW 
cutoff. For diamond, for example, the combination of 



3.4 Oq and G 



4.2 Oq leads to stable re- 



C 

^max 

suits. 

An overall view of the MT and the interstitial ex- 



change potential on the diamond (Oil) plane is shown 
in Fig. mja) as a contour plot. The plane is displayed in 
Fig. [21 It contains the connecting lines along [111] and 
JlOO] that correspond to Fig.[T] We see that in the regions 
close to the atomic nuclei the potential is predominantly 
spherical. However, towards the MT sphere boundaries 
the potential becomes strongly anisotropic and matches 
continuously to the warped interstitial potential, which 
is far from constant also. In fact, the nonsphericity of the 
EXX potential is considerably more pronounced than in 
the LDA potential, Fig. Hfb). The latter is similar in 
shape to the electron density distribution Jcf. Fig. [5ja)], 
of which it is a direct function Kc^°^(r) oc nirf/^. The 
EXX potential, in contrast, incorporates the full nonlo- 
cality of the EXX functional, which makes it much more 
corrugated than the LDA one, in particular in the MT 
spheres, where the KS orbitals are highly oscillatory. All 
this stresses the importance of a full-potential treatment 
within the EXX-OEP approach. 

The LDA potential corresponds in each point r to the 
exchange potential of the homogeneous electron gas with 
an electron density that equals the local electron density 
n(r) of the real system. Thus, by construction it is exact 
for the homogeneous electron gas but misses the effects of 
density variations. The EXX potential, in contrast, takes 
all density variations exactly into account. Thus, the 
differences between Figs.Hl^a) and (b) must be attributed 
to the influence of the density inhomogeneities on the 
exchange potential. This influence is particularly large 
in regions where the density varies a lot, that is, close to 
the atomic nuclei, while in the interstitial region the two 
potentials are more similar. 

The differences in the exchange potentials naturally af- 
fect the electron density distribution. The EXX electron 
density. Fig. Efa), clearly shows a pronounced contrac- 
tion of the electron distribution compared with the LDA 
one. This is a direct consequence of the self-interaction 
error, which is eliminated in the EXX approach, while 
it gives rise to an unphysical delocalization in the case 
of the LDA potential. This becomes clearer in Fig. [Sfb) 
where we plot the difference between the EXX and the 
LDA densities. The exactly compensated self-interaction 
allows the charge to accumulate in the atomic cores, but 
also in the covalent bonds between the atoms. 

To understand the requirement of the basis-set balance 
in more detail, we go back to the OEP Eq. ([T8l) . whose 
solution involves the inversion of the response function, 
Eq. The response function describes the linear re- 

sponse of the electron density with respect to changes 
of the effective potential. For the latter we employ the 
MPB, while the former is given by the orbital densities 
of the occupied states, Eq. (|4]), and, hence, ultimately 
by the LAPW basis set. Thus, the LAPW basis must 
provide enough flexibility for the density to enable it to 
respond adequately to the changes of the effective poten- 
tial. 

This explains the observed behavior and becomes evi- 
dent in the convergence of the response function with re- 



(a) 



(a) 




(b) 

Figure 4: Two-dimensional plot of (a) the local EXX and (b) 
the LDA exchange potential on the (Oil) plane of diamond. 
Both potentials are shifted such that fV^{r)d^r = 0. The 
contour lines start at 14 = — l.OOhtr and have an interval of 
0.05 htr. The contour line with the maximal value is found 
(a) at K = 0.30 htr and (b) at K = 0.20 htr. The dotted line 
corresponds to 14 = 0. The MT sphere boundaries and the 
lines corresponding to Fig. [1] are indicated. 



spect to the LAPW basis. In Fig.lHlwe show the changes 
of the eigenvalues of Xsi ordered according to increas- 
ing moduli, when we add more and more local orbitals, 
21 + 1 additional local orbitals per / quantum number 
{I = 0, . . . , 5) in each step. We again use L,„ax = 4 and 




(b) 

Figure 5: (a) Total electron densities obtained from the 
EXX (solid lines) and LDA potential (dotted lines) on 
the diamond (Oil) plane. Contour lines for the values 
0.03 flo^ 0.06 tto^..., 0.99 flo^ are shown, (b) Plot of the 
density difference with contour lines between —0.039 a^^ and 
0.390 a,^'^ at intervals of 0.003 a^^. The dotted line corre- 
sponds to An(r) = 0. 



G'max — 3.4 Oq^ as MPB cutoff values. In the case of the 
maximal number of local orbitals per Im channel, nio = 6, 
the basis is increased by 432 functions relative to the con- 
ventional LAPW basis. Clearly, the eigenvalues can be 
systematically converged. The relative changes are in the 
order of 0.1% to 1.0% between n\o = 5 and njo — 6. Es- 
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Figure 6: Convergence of the eigenvalues Xs,i of the response 
function for diamond. The Xs,i are all negative and or- 
dered according to increasing moduli. The relative change 
[Xs,-f ('^lo + 1) -Xs,-r('^lo)]/Xs,/('^lo) is plotted, where nio is the 
number of local orbitals per Im channel {I — 0, ... ,5, \m\ < I). 
The case nio — corresponds to the conventional LAPW ba- 
sis without local orbitals. A Bezier algorithm was employed 
to smooth the curves. 



pecially the small eigenvalues converge well, which are 
particularly important in the inversion of the response 
function. We note that a straightforward elimination of 
the small eigenvalues by singular value decomposition is 
not advisable and leads to an ill-defined response func- 
tion. 

Up to now, we have discussed the importance of the 
quality of the LAPW basis for the shape of the local EXX 
potential. Only for well-balanced basis sets, the poten- 
tial is smooth and physical (cf. Fig. [T|). We now address 
the question to what extent this effect influences the KS 
one-particle energies that result from the self-consistent 
solution of Eq. ^ with the local EXX potential obtained 
from the OEP Eq. Table HVl gives the transition 

energies, that is, the KS eigenvalue differences, from the 
valence-band maximum at the T point to the L and X 
point of the lowest conduction band of diamond for the 
basis sets with nio = 0, . . . , 6 local orbitals per Im chan- 
nel. Obviously, at least three local orbitals are necessary 
to converge these transition energies to within 0.01 eV. 
Between the unbalanced and balanced basis sets (nio = 
and nio — 6, respectively) the values change by about 0.2 
eV. 

Another balance condition we find for the Is core state 
that goes into both the left- and the right-hand sides of 
Eq. ((T8| . In particular, we now distinguish between four 
cases: (a) core state considered in Xs,/j and tj, (b) core 
state only considered in Xs,/jj (c) core state only consid- 
ered in t], and (d) core state neglected in both. Case 
(a) corresponds to the full calculations presented so far. 
Figure [7] shows that the resulting potentials look very dif- 
ferent for the different cases. Surprisingly, potential (d) is 
closest to the full potential, while the potentials (b) and 





r ^ r 


r ^ L 


r X 





6.351 


9.243 


5.307 


1 


6.196 


9.086 


5.125 


2 


6.186 


9.069 


5.144 


3 


6.180 


9.063 


5.138 


4 


6.178 


9.059 


5.139 


5 


6.177 


9.057 


5.136 


6 


6.176 


9.055 


5.136 



Table IV: KS transition energies (in eV) for diamond obtained 
from the self-consistent solution of the KS equation with the 
local EXX potential for LAPW basis sets including zero to 
six local orbitals per Im channel (Z = 0, . . . , 5, |m| < I). As 
the LAPW basis is made more flexible, the transition energies 
converge. 




0.5 1 1.5 2 2.5 



Figure 7: Same as Fig. [T] for the cases (a) Is state fully con- 
sidered, (b) only considered in Xs,iJ, (c) only considered in 
ti, and (d) neglected in both (see text). 



(c) are much too shallow and too strongly varying, re- 
spectively. Obviously, the inclusion of the core state only 
on one side of the OEP equation gives rise to an equa- 
tion that is out of balance and that yields an unphysical 
V^(r). This also influences the resulting KS transition 
energies. In the balanced case (d), the results 6.19 eV, 
9.08 eV, and 5.28 eV for F ^ F, F ^ L, and F -> X, 
respectively, are surprisingly close to the full calculation 
(a) (cf . Table lIVp , while the energies for the unbalanced 
case (b) deviate more strongly, especially for the F X 
transition, 6.22 eV, 9.16eV, and 6.03 eV. The energetic 
position of the Is state with respect to the Fermi en- 
ergy, however, is only realistic for the full calculation (a), 
e(ls) ^ -265.7 eV (cf. -262.9 eV for LDA), whereas (b) 
and (d) give binding energies, that are smaller by 31.6 eV 
and 10.5 eV, respectively. Calculation (c) is unstable and 
does not converge. These results indicate that the PP 
approximation is, indeed, suitable for the EXX-OEP ap- 
proach, a conjecture that will be confirmed by our refer- 
ence calculations later-on. 

As outlined in Sec. |lTl the construction of a local EXX 
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potential within the KS formahsm of DFT is equivalent 
to the OEP approach of Sharp and Horton^i^ where the 
Hartree-Fock total energy is minimized under the con- 
straint that the wave functions feel a local multiplicative 
potential. This constraint reduces the Hilbert space for 
the wave functions and, thus, increases the total energy 
due to the variational principle. In fact, we find that 
the total energy for diamond obtained with the nonlocal 
Hartree-Fock potential, Eq. ([TU)) . and an 8x8x8 k-point 
set is 0.24 eV per unit cell lower than that of the EXX- 
OEP approach. 

As reference, we now report fully converged KS transi- 
tion energies for a variety of semiconductors and insula- 
tors in Table IVl All calculations are performed with an 
8x8x8 Brillouin zone sampling and at the experimental 
lattice constants (C: 6.746 cq, Si: 10.26 ao, SiC: 8.24 ao, 
Ge: 10.67 ao, GaAs: 10.68 oq, Ne: 8.44 ao, Ar: 9.93 ao). 
Apart from the EXX-only calculations, we also show 
transition energies obtained with the EXX+VWNc func- 
tional in which the LDA correlation functional from 
Ref. 1^ was added to the EXX functional. Table |V] shows 
that both the EXX and EXX+VWNc functionals yield 
KS transition energies much closer to experiment (last 
column) than the LDA functional (first column). While 
for semiconductors the resulting energies are even quan- 
titatively in very good agreement with experiment, there 
are larger discrepancies for the insulators diamond and, 
in particular, crystalline Ne and Ar. There is only lit- 
tle difference between the EXX and EXX+VWNc val- 
ues. The inclusion of the LDA correlation functional 
does not lead to a definite improvement. The direct 
band gap F — T for diamond of 6.21 eV agrees very 
well with the value 6.18 eV recently reported by Engel,^^ 
who used a plane-wave PP approach with cutoff values 
pushed to the AE limit. It is even identical to the value 
obtained with a standard valence-only plane-wave PP 
approach.^- In contrast to that, Sharma et al. calcu- 
lated a much larger value of 6.67 eV with their FLAPW- 
EXX implementation. For neon, there is a somewhat 
larger discrepancy with the calculation by Magyar et 
al.^ though. In conclusion, with our EXX-OEP im- 
plementation within the AE FLAPW method, we obtain 
results in very good agreement with previous plane- wave 
PP calculations, provided that the basis sets for the wave 
functions and the potential are properly balanced. This 
shows that the PP approximation is adequate for the 
EXX-OEP approach at least for the systems examined 
here, which is at variance with the findings of Ref. 25. 



VI. CONCLUSIONS 

We have developed an all-electron full-potential im- 
plementation of the EXX-OEP approach to DFT within 
the FLAPW method. We analysed the conditions and 
requirements on the basis sets and numerical cutoff pa- 
rameters to obtain reliable and numerically stable results. 
Based on this knowledge we presented as proof of prin- 



ciple results on KS transition energies for some typical 
semiconductors, insulators and noble-gas solids that are 
in very good agreement with pseudopotential results. 

The OEP equation is formulated utilizing the mixed 
product basis (MPB),"*^"^"* which has been adjusted for 
the present purpose: it is augmented with the atomic 
EXX potential, the constant function is eliminated, and 
the basis functions are made continuous all over the 
space. In this basis, the OEP equation becomes an alge- 
braic equation, which is solved for the local EXX poten- 
tial with standard numerical tools. 

For the case of diamond, we have demonstrated that 
the local EXX potentials are spatially strongly corru- 
gated, which makes a full-potential treatment even more 
important than in conventional LDA or GGA calcula- 
tions. Furthermore, the two basis sets, LAPW and MPB 
are not independent. They must be properly balanced 
to obtain a smooth and physical EXX potential over the 
whole space. In the unbalanced case, the potential shows 
spurious oscillations, which we have traced back to an in- 
sufficiently converged response function, a function that 
gives the response of the electron density with respect to 
changes of the effective potential. If the LAPW basis, 
which parametrizes the electron density, is not flexible 
enough, the electron density cannot follow the changes 
of the effective potential that are described by the MPB 
leading to a corrupted response function. As a result, 
the LAPW basis must be converged with respect to a 
given MPB. Already in the simple case of diamond, we 
must add six local orbitals at different energies for each 
Im channel from ^ = to ^ = 5 in order to obtain a 
smooth potential in the spheres. This shows that the 
LAPW basis must be converged to an accuracy that is 
far beyond that of conventional LDA or GGA calcula- 
tions. Similarly, also the LAPW reciprocal cutoff radii 
must be chosen large enough. 

Not surprisingly, the shape of the EXX potential - os- 
cillatory or smooth - has an impact on the resulting KS 
transition energies. We find that with properly balanced 
basis sets, the transition energies for a variety of semi- 
conductors and insulators obtained with the EXX and 
the EXX+VWNc functionals are in very good agree- 
ment with plane-wave pseudopotential results from the 
literature (crystalline neon is an exception). This con- 
firms that the pseudopotential approximation works re- 
liably within the EXX-OEP approach. Our finding is in 
contradiction to a previously published implementation 
(Ref. [25I ) based on the FLAPW method, where large dis- 
crepancies with pseudopotential results were reported. 

Currently, reliable all-electron full-potential EXX- 
OEP calculations are computationally very demanding, 
because of the need for large orbital basis sets, which we 
attribute partly to the fact that the LAPW basis func- 
tions depend explicitly on the effective potential. To re- 
fine our full-potential implementation of the EXX-OEP 
approach, we suggest as an task for the future the inves- 
tigation of schemes that treat the response of the LAPW 
basis with respect to changes of the potential more effi- 
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Table V: KS transition energies (in eV) obtained with the local EXX and EXX+VWNc potentials and an 8x8x8 k-point 
sampling. For comparison, plane-wave PP results and experimental values from the literature are given. 



ciently than employing local orbitals. 
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